(* Mathematica package *)

(*Internally implemented definitions of SI and SIIv, less pattern matching required but no argument error checks.  These
provide results for subsystems up to when the are effectively S so that Phiv and the alternate partial system transfer matrix are
taken as 2 x 2 identity matrices*)

AbelesSIv[pol_, omega_, thetalist_, nlist_, Dlist_, v_] /; v == 0 := IdentityMatrix[2];
AbelesSIv[pol_, omega_, thetalist_, nlist_, Dlist_, v_] /; v == 1 := 
 RefractionMv[pol, thetalist[[1]], thetalist[[2]], nlist[[1]], nlist[[2]]];
 
AbelesSIv[pol_, omega_, thetalist_, nlist_, Dlist_, v_] /; 2 <= v <= Length@Dlist+1 := 
	Dot @@ Riffle[
    RefractionMv[pol, thetalist[[1 ;; v]], thetalist[[2 ;; v + 1]], nlist[[1 ;; v]], nlist[[2 ;; v + 1]]], 
    PhaseMv[omega, thetalist[[2 ;; v]], nlist[[2 ;; v]], Dlist[[1 ;; v - 1]]]
    ];

(*More restrictive pattern-match for args, for use by users, executes just about as fast as internal definitions. These
execute ~ 1.1 - 1.5x longer than internally implemented SIIv*) 
AbelesSIv[pol_?PolQ, omega_, ThetavNvDv : {{_, _}, {_, _, _} .., {_, _}}, v_Integer] /; v > Length@ThetavNvDv - 1 :=
	AbelesS[pol, omega, ThetavNvDv]

AbelesSIv[pol_?PolQ, omega_, ThetavNvDv : {{_, _}, {_, _, _} .., {_, _}}, v_Integer] /; 2 <= v <= Length@ThetavNvDv - 1 :=
 Dot @@ Riffle[
   RefractionMv[pol, ThetavNvDv[[1 ;; v, 1]], ThetavNvDv[[2 ;; v + 1, 1]], ThetavNvDv[[1 ;; v, 2]], ThetavNvDv[[2 ;; v + 1, 2]]], 
   PhaseMv[omega, Sequence @@ Transpose[ThetavNvDv[[2 ;; v]]]]
   ];
AbelesSIv[pol_?PolQ, omega_, ThetavNvDv : {{_, _}, {_, _, _} .., {_, _}}, v_Integer] /; v == 1 :=
  RefractionMv[pol, Sequence @@ Flatten@Transpose[ThetavNvDv[[1 ;; 2, 1 ;; 2]]]];
AbelesSIv[pol_?PolQ, omega_, ThetavNvDv : {{_, _}, {_, _, _} .., {_, _}}, v_Integer] /; v == 0 := IdentityMatrix[2];


AbelesSIv[pol_?PolQ, omega_, Theta0_, n0_, NvDv:{{_,_}..}, nkp1_, v_Integer] /; v > Length@NvDv :=
	AbelesS[pol, omega, Theta0, n0, NvDv, nkp1];
AbelesSIv[pol_?PolQ, omega_, Theta0_, n0_, NvDv:{{_,_}..}, nkp1_, v_Integer] /; 2 <= v <= Length@NvDv :=
	Block[{ThetaLIST = {Theta0} ~Join~ RefractionAngle[Theta0, n0, NvDv[[All,1]]~Join~{nkp1}],
	           nLIST = {n0} ~Join~ NvDv[[All,1]] ~Join~ {nkp1}},
	   Dot @@ Riffle[
		RefractionMv[pol, ThetaLIST[[1 ;; v]], ThetaLIST[[2 ;; v + 1]], nLIST[[1 ;; v]], nLIST[[2 ;; v + 1]]], 
		PhaseMv[omega, ThetaLIST[[2 ;; v]], nLIST[[2 ;; v]], NvDv[[1 ;; v - 1, 2]]]]];
AbelesSIv[pol_?PolQ, omega_, Theta0_, n0_, NvDv:{{_,_}..}, nkp1_, v_Integer] /; v == 1 :=
  RefractionMv[pol, Theta0, RefractionAngle[Theta0, n0, NvDv[[1,1]]], n0, NvDv[[1,1]]];
AbelesSIv[pol_?PolQ, omega_, Theta0_, n0_, NvDv:{{_,_}..}, nkp1_, v_Integer] /; v == 0 := IdentityMatrix[2];

AbelesSIv[___]/;Message[General::badargs,AbelesSIv]:="dummy result, so symbol returned unevaluated";

